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Abstract. Recent experiments on cortical neural networks have revealed the 
existence of well-defined avalanches of electrical activity. Such avalanches have been 
claimed to be generically scale-invariant - i.e. power-law distributed - with many 
exciting implications in Neuroscience. Recently, a self-organized model has been 
proposed by Levina, Herrmann and Geisel to justify such an empirical finding. Given 
that (i) neural dynamics is dissipative and (ii) there is a loading mechanism "charging" 
progressively the background synaptic strength, this model/dynamics is very similar 
in spirit to forest-fire and earthquake models, archetypical examples of non-conserving 
self-organization, which have been recently shown to lack true criticality. Here we 
show that cortical neural networks obeying (i) and (ii) are not generically critical; 
unless parameters are fine tuned, their dynamics is either sub- or super-critical, even if 
the pseudo-critical region is relatively broad. This conclusion seems to be in agreement 
with the most recent experimental observations. The main implication of our work is 
that, if future experimental research on cortical networks were to support that truly 
critical avalanches are the norm and not the exception, then one should look for more 
elaborate (adaptive/evolutionary) explanations, beyond simple self-organization, to 
account for this. 
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1. Introduction and outlook 

1.1. Generic scale invariance 

In contrast to what occurs for standard criticality, where a control parameter needs to 
be carefully tuned to observe scale invariance, certain phenomena as earthquakes, solar 
flares, avalanches of vortices in type II superconductors, or rainfall, to name but a few, 
exhibit generic power-laws, - i.e. they lie generically at a critical point without any 
apparent need for parameter fine tuning [HE]. Ever since the concept of self-organized 
criticality [1] was proposed to account for phenomena like these, it has generated a lot 
of excitement, and countless applications to almost every possible field of research have 
been developed. Underpinning the necessary and sufficient conditions for a given system 
to self-organize to a critical point is still a key challenge. 

In this context, it has been established from a general viewpoint that conserving 
dynamics (i.e. that in which some quantity is conserved along the system evolution) is 
a crucial ingredient to generate true self-organized criticality in slowly driven systems 
[31 H]. In this way, non-conserving self-organized systems have been shown not to be 
truly scale- invariant (see [1] and references therein). While sandpiles, ricepiles, and other 
prototypical self-organized models are examples of conserving self-organizing systems, 
forest-fire and earthquake automata are two examples of non-conserving models. They 
both were claimed historically to self-organize to a critical point and they both were 
shown afterwards to lack true scale- invariant behavior (see [1] and references therein). 
The main reason for this is, in a nutshell, that non-conserving systems combine driving 
(loading) and dissipation, and this suffices to keep the system "hovering around" a 
critical point separating an active from a quiescent or absorbing phase (driving slowly 
pushes the system into the active phase and dissipative takes it back to the absorbing 
phase). But, in order to have the system lying exactly at the critical point requires of an 
exact cancellation between dissipation and driving (loading) to achieve a critical steady 
state; such a perfect balance can only be achieved by parameter fine tuning, and then 
the system cannot be properly called "self-organized" . 

This mechanism of (non-conserving) self-organization has been termed self- 
organized quasi-criticality (SOqC) [4] to underline the conceptual differences with truly 
scale-invariant, (conserving) self-organized criticality (SOC) [H [2]. From now on, we 
shall use the acronym SOqC to refer to non- conserving self-organized systems, and shall 
keep the term SOC for self-organized conserved systems. 

SOqC may explain the "approximate scale invariance" (with apparent power-law 
behavior extending for a few decades) observed in many real systems as those mentioned 
above (earthquakes and forest fires) but, strictly speaking, it fails to explain true scale- 
invariance. SOqC systems require some degree of parameter tuning to lie sufficiently 
close to criticality. For a much more detailed explanation of the SOqC mechanism and 
its differences with SOC, we refer the reader to [1]. 
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1.2. Scale invariance in neuronal avalanches? 

Neuronal avalanches were first reported by Beggs and Plenz, who analyzed in vitro 
cortical neural networks using slices of rat cortex as well as cultured networks O El [7] . 
More recently, neuronal avalanches have been observed also in vivo [8J. In all these 
cases, cortical neurons form dense networks which, under adequate conditions, are able 
to spontaneously generate electrical activity [9J. The associated local field potentials can 
be recorded by using multielectrode arrays |5J. Each electrode in the array monitors 
the electrical activity of a local group of neurons (which for convenience can be thought 
of as a unique "effective" neuron; a review of the involved experimental techniques 
and methods can be found in [I^). According to Beggs and Plenz [5l El [7j activity 
appears in the form of "avalanches" , i.e. localized activity is generated spontaneously 
at some electrode and propagates to neighboring ones in a cascade process which 
occurs at a much faster timescale (tens of milliseconds) than that of the quiescent 
periods between avalanches (typically of the order of seconds). Previous experimental 
research in cultured networks had identified the existence of spontaneously generated 
synchronized bursts of activity (involving synchronous activation of many neurons), 
followed by silent periods of variable duration [HI [121 1131 HH US] (theoretical work 
has been done to explain such a coherent or synchronous behavior, see for instance, 
[T6l [T7]). The main breakthrough by Beggs and Plenz in [5jj was to enhance the 
resolution and bring the internal structure of "synchronized" bursting events to light. 
In other words, the apparently synchronous activation of many neurons required for 
a synchronized burst corresponds to a sequence of neuron activations, i.e., a neuronal 
avalanche, which generates spatio-temporal patterns of activation confined between two 
consecutive periods of quiescence. 

Experimental measurements of avalanches can be performed, and the distribution 
of quantities as i) the avalanche size s (i.e. the number of electrodes at which a non- 
vanishing signal is detected during an avalanche) and ii) the avalanche lifetime, t, can be 
recorded. What is relevant for us here is that, according to Beggs and Plenz, avalanches 
seem to be generically scale invariant [51E1IZ]; in particular, avalanche sizes, s, and times 
t are distributed as: 

p{s) ~ s-'/'T{s/s,), p{t) ~ t-'g{t/t,) (1) 

respectively, where J-" and Q are two cut-off functions; the cutoff Sc grows in a scale 
invariant way as a function of system-size: the larger the system the larger the cutoff, 
providing evidence for finite size scaling. The cut-off tc appears at very small times, so 
the evidence for scale invariance is much larger for s than for t. 

These results have been claimed to be robust across days, samples, and 
pharmacological variations of the culture medium [HI El E]- The exponent values in 
Eq. ([T|) coincide with their mean-field counterparts for avalanches in sandpiles (the 
prototypical examples of self-organized criticality) [18]. Mean- field exponents do not 
come as a surprise: given the highly entangled structure of the underlying network 
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(which has been reported to have the small-world property [IS]) mean-field behavior is 
to be expected for critical phenomena occurring on it [20]. 

Finally, recalling that, at a mean-field level, avalanche dynamics can be interpreted 
as a branching process [21] , an empirical study of the branching ratio, a, (defined as the 
fraction of active electrodes per active electrode at the previous time bin) was performed 
in [5l |22]. It was found that the value of a measured for avalanches started from one 
single electrode is very close to unity, in agreement with the critical value of marginally 
propagating branching processes, dc = 1. 

From these results, it has been claimed that cortical neural networks are generically 
critical, i.e. scale-invariant, and that they reach such a critical state in a "self-organized" 
way [5]. Scale invariance in the propagation of neural activity has raised a great deal of 
interest and excitement in Neuroscience. For instance, critical neural avalanches have 
been claim to lead to [I2l [5] : 

• optimal transmission and storage of information [5l [6l [TJ |22l [23] , 

• optimal computational capabilities [21], 

• large network stability [25] , 

• maximal sensitivity to sensory stimuli [26], etc. 

Let us caution that discrepant results, i.e. non-critical neuronal avalanches, have 
also been recently reported in the literature. For instance, measurements of cortical 
local-field-potentials were performed by Bedard et al. [27] using parietal cat cortex. 
None of the features reported by Beggs and Plenz [5] was observed for such a network; 
not only the observed behavior was not critical, but it was not even possible to observe 
clean-cut avalanches. It was argued that the absence of scale-free avalanches could stem 
from fundamental differences between the considered cortex regions used in [27] and in 
[5] . Moreover, in a recent review paper, Pasquale et al. [28] report on different empirical 
types of avalanche distributions: critical, subcritical, or super-critical, depending on 
various factors. These authors conclude that critical avalanches can indeed emerge, but 
they are more the exception that the rule. 

1.3. Goals and outlook 

The main goal of this paper is to elucidate from a theoretical viewpoint whether neuronal 
avalanches are truly critical or not. Or, more precisely, to understand whether self- 
organizing mechanisms (such as those of SOC or SOqC) can justify the findings for 
neuronal avalanches. To this purpose, we rely extensively on a model for neuronal 
avalanches, proposed recently by Levina, Herrmann and Geisel [29] . The model is a self- 
organized one, including integrate-and-fire neurons and short-term synaptic plasticity. It 
has been claimed, both analytically and numerically, to back the existence of generically 
(strictly) critical neuronal avalanches in a very broad region of parameter space [29] . 

The key observation, which motivated the present work, is the fact that local 
conservation laws, such as those required to have truly critical self-organized (SOC) 
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behavior, are not present in neural networks in any obvious way. If cortical networks 
are represented as an electrical circuit, perfect transmission without loss of energy is an 
unrealistic idealization and, analogously, if they are modeled as networks of dynamical 
synapses, there also exist dissipative or "leakage" phenomena. In summary, no quantity 
is strictly conserved in neural signal transmission. Reasonably enough, the Levina, 
Herrmann and Geisel (LHG) model [29] is also a non-conserving one (see below). 

Therefore, the existence of critical neuronal avalanches (both experimentally and 
in the LHG model) seems to be in contradiction with the general conclusion in [1], i.e. 
the lack of true criticality in non-conserving systems. In this way, a rationalization of 
neuronal avalanches would only be possible, at most, in terms of self-organized quasi- 
criticality (SOqC) and not in terms of strict criticality as suggested in [29] . 

Following the steps in [4], here we shall underline the analogies and differences 
between the model by Levina et al. and other non-conserving self-organized models 
such as those for earthquakes or forest fires. We shall show that the LHG model is 
not generically critical: it can be either critical, subcritical or supercritical depending 
on parameter values; fine tuning is required to achieve strict scale- invar iance. Still, the 
model is capable of generating, for a relatively wide parameter range, pseudo-critical 
avalanches with associated truncated power-laws which can suffice to explain empirical 
observations. 

This conclusion, -i.e. the lack of true criticahty- is expected to apply not only to 
the model in [29], but also to empirical neuronal avalanches. It suggests that if neuronal 
avalanches turned out to be truly critical, the ultimate reason for that should be looked 
for in some type of adaptive/evolutionary mechanism [30] or in homeostatic processes 
[3T] . but cannot be generically ascribed to plain self-organization. 

The rest of the paper is structured as follows. In Section [2l we present the self- 
organized model proposed by Levina et al. for neuronal avalanches. A discussion of its 
main properties appears in Sections |3] (numerical) andH] (analytical). Then, in Section 
O we put this model under the general framework of self-organized quasi-criticality 
introduced in [1] by deriving explicitly a Langevin equation from its microscopic rules, 
emphasizing the lack of true generic criticality. Finally, the main conclusions and a 
critical discussion of recent experimental results are presented. 

2. The Levina-Herrmann- Geisel (LHG) model 

Aimed at understanding the origin of power-law distributed cortical avalanches, Levina, 
Herrmann, and Geisel (LHG) [29] proposed a variation of the well-known Markram- 
Tsodyks model of chemical synapses [17]. Such a model had been already extensively 
used to reproduce the dynamics of synchronized bursting events (also called "population 
spikes") [13132]. 

Consider a fully connected network of N integrate-and-fire neurons each of them 
characterized by its local (membrane) potential, Vi, with 

0<V,< Vmax. (2) 
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Neurons i and j (with i ^ j) are connected by a synapse of strength Jij. This can be 
thought of as the amount of available neurotransmitters or, more generally, "synaptic 
resources" , for such a connection. 

In the original Markram-Tsodyks model [17], together with Vi and Jij, there is 
a third variable, Uij, representing the fraction of neurotransmitters which is actually 
released every time a pulse is transmitted between i and j. Its dynamics can be used to 
implement synaptic facilitation (see, for instance, [33]); but, aimed at keeping the model 
as simple as possible, and following LHG [29], we fix Uij = m to be a constant. 

The simplified Markram-Tsodyks or LHG dynamics is defined by the following 
equations: 

-i = j-*5(t - t^j + J2w^^(^ - ^sp) - y^aM - tl^) 

3=1 ~ ^ (3) 



dJij 1 /a 



dt Tj \U 'V ' 

The different terms in Eq. follows: 

• Driving: J^^'* is the amplitude of an external random input which operates at 
discrete times f^^^^ on i. Driving impulses can be introduced at a fixed rate h. 
Alternatively, slow-driving {h — t- oo) can be implemented by switching on if 
and only if all potentials are below threshold. 

• Firing: —Vmax^it — tip)] if the potential at i overcomes the threshold, Vmax, at time 
tip, the neuron spikes, and it is reset to 

V^itlp) ^ VMp) - V^ax; (4) 



otherwise, nothing happens. 

N-l 

• Integration: ^^""^'^(^ ~ ^sp)) the (post-synaptic) neuron i integrates signals of 

j=i ~ 

amplitude uJij/ {N— 1) from each spiking (pre-synaptic) neuron j. A non- vanishing 
delay between the time of discharge and the time of integration in neighboring 
neurons could also be introduced, without affecting significantly the results. 

• Synaptic depression: —uJij6(t — tip); after each discharge involving the (pre- 
synaptic) neuron j all synaptic strengths Jij (where i runs over all post-synaptic 
neurons) diminish by a fraction u. 

• Synaptic recovery: — ( Jj , V synapses recover to some target value, Jij = J = 

Tj \u / 

on a timescale determined by the recovery time, tj. 

Observe that the only sources of stochasticity are the initial condition and the 
external driving process, while the avalanche dynamics is purely deterministic. Also, 
the set of equations above can be implemented on any generic network topology; here 
(following LHG) we will mostly restrict ourselves to fully connected networks, even if 
results for random networks and two-dimensional lattices are also briefly discussed. 
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3. Model analysis 

3.1. Static limit 

Let us first discuss the static limit of tlie model in which the synaptic recovery rate is so 
fast (i.e Tj — )■ 0) that Jij can be taken as a constant for all pairs i,j and for all times: 
•^hj = J = a''*"**'^/^. In such a case (keeping u fixed), ct'^*"**'^ acts as a control parameter 
[34] . Observe that, in the limit in which q/**'^*'^ — ). Vmaxi the model becomes conserving: 
each spiking neuron reduces its potential by Vmax and each of its {N — 1) neighbors is 
increased by Vmax/{.N — 1) (integration term in Eq. (j3])). 

Once the system has reached its steady state, it is possible to assume that the 
values of V are uniformly distributed in the interval [e, Vmax — e] with e — )■ when 
— 7- oo. This assumption parallels what is done in a similar analysis of related self- 
organized systems such as earthquake models |35] and can be numerically verified to 
hold with good accuracy (see Appendix). This implies that, fixing (without loss of 
generality) Vmax = 1, in the large system-size limit, a randomly chosen neuron can be 
in any possible state with uniform probability. Thus, upon receiving a discharge of size 
uJ / {N — 1), it becomes over threshold with probability uJ/ {N — 1). Hence, viewing the 
propagation of activity within avalanches as a branching process with branching rate 
uJ/{N — 1) and — 1 neighbors per neuron, the average avalanche size (s) can be 
written as the sum of an infinite geometric series [2T] 



Observe that this expression is valid only for uJ < 1. The model critical point can be 
identified by the presence of a divergence in Eq. ([5]); this occurs at the conserving limit 
^static _ agreement with what happens in other models of SOC (like sandpiles) 
which are critical only in the case of conserving dynamics. 

For a**"**'^ > 1 (i.e. above the conserving limit) the potential at each site grows 
unboundedly (i.e. there is no stationary state) with perennial activity (generating an 
"explosive" super-critical phase) while, for q/''*'***'= < 1^ the process is dissipative on 
average, i.e. the total potential is reduced at every spike and avalanches die after a 
characteristic time (sub-critical phase). Thus, in summary, as already discussed in the 
literature [3l], the static version of the LHG model exhibits a standard (absorbing) phase 
transition separating a sub-critical from a super- critical phase. 

Let us remark that, for finite systems, the critical point has size- dependent 
corrections. It is only in the infinite size limit, in which driving and dissipation 
vanish, that q/^*'***'= = uJc = 1. Actually, for any finite system, e 7^ 0, and additional 
finite-size terms need to be included in the calculation above. This is a consequence 
of the fact that, in order to achieve a steady state for finite systems, some form of 
dissipation needs to be present to compensate the non- vanishing driving, J^^*, entailing 
^siaticj^jy-j < a^*°**'^(A^ — )• oo) = 1 (scc Table 1 where numerical estimates for the critical 
point location are reported; details of the computational procedure are reported in the 




(5) 
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2000 
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^static 


0.92(1) 


0.93(1) 


0.94(1) 


0.95(1) 


0.96(1) 


0.97(1) 




1 



Table 1. Location of the critical point a^'"**'^ as a function of the system size N, 
as obtained in computer simulations of the static model (tj — >■ 0). The critical point 
location does not depend on the way the system is driven, i.e. on /^^*. 

forthcoming section). 
3.2. Dynamic model 

Let us now turn back to the full dynamic model. Observe that: 

• The equation for Jjj in Eq. (|3]) includes a loading mechanism (analogous to 
those reported in [1] for earthquake and forest-fire models) or "synaptic recovery 
mechanism" which counterbalances the effect of synaptic depression in the absence 
of spikes: the "background field" Jij increases steadily towards its target value 
a/u. Note that in contrast with models of forest fires or earthquake automata, in 
which the "loading mechanism" (see [1]) acts only between avalanches, the recovery 
djTiamics of J occurs also during avalanches, at a finite timescale controlled by tj. 

• In the limit in which u{Jij) — )■ Vmax "^{hj), where (.) stands for steady-state time 
averages, conservation is recovered on average. In analogy with the static model, the 
dynamics becomes non-stationary above such a limit: loading overcomes dissipation 
and potential fields grow unboundedly. 

In the case in which lext drives the system slowly we are in the presence of the 
main ingredients characteristic of non-conserving self-organized models, as described 
generically in |1]: 

(i) separation of (driving and dynamics) timescales, 

(ii) dissipative dynamics (provided that u{Jij) < Vmax), and 

(iii) loading mechanism, increasing the average value of the "background field" Jij. 

Prior to delving into further analytical calculations, which are left for Section Hj let 
us present in the rest of this section computational results obtained for Eq. (|3]). 

3.2.1. Numerical analyses Numerical integration of Eq. (|3]) becomes very costly as 
the number of components grows, limiting the maximum system size (up to = 3000 
in the present study). Observe that, owing to the presence of 5-functions, Eq. ([3]) is 
an "impulsive dynamics" equation and thus, caution must be paid when integrating it 
numerically not to miss delta peaks when discretizing. 

The system is initialized with arbitrary (random) values of G [0, Vmax] and 
Jjj G [0,1] W{i,j). We keep a as a control parameter and fix parameter values 
mostly as in [29]: u = 0.2, Vmax = 1 and tj = ION. 
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Figure 1. Time evolution of uJ and tlie number of spiking neurons for a — 0.9 
subcritical (up), a — lA critical (center), and a = 1.9 supercritical (down), in 
simulations with N = 1000. It is only above the critical point of the dynamical model 
that uJ goes beyond the critical point of the static model for the considered system 



size, af"*^" = 0.95(1) (dashed line). 



Let us remark that the choice tj = ION [29] might be not very reahstic from a 
neuro-scientific point of view; i.e. it is not clear whether the synaptic recovery rate 
should depend on the total number of connections of the corresponding neuron or not. 
Observe that N is the number of synapses per neuron, therefore in principle, it could 
be the case that, if a given neuron has limited resources, the recovery rate per synapse 
depends on the total number of synapses. But also the opposite could be true; i.e. the 
recovery time of a given synapse could depend only on its local properties and not on 
those of its corresponding neuron. This is a neuro-scientific issue that is beyond the 
scope of the present paper and that we prefer not to enter here. Anyhow, we have 
verified that our results are not significantly affected by such a choice; for example, 
we have also considered values of r fixed for any and checked the robustness of our 
results. 

We work in the slow driving limit, i.e. we drive the system with an input, lext at 
a randomly chosen site if and only if all potentials are below threshold. The sequence 
of activity generated therefrom constitutes an avalanche. We have used two different 
A^-dependences for /^^*: (i) /^^* = 7.5 x A^-^, and (ii) J^^'* = A^-o-6467. ^q^j^ ^-^^^ 
engineered to comply with the scaling form f^^'* ~ A^""" considered by Levina et al. , and 
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Figure 2. Left: Avalanche-size distribution of the LHG model for A'^ = 1000 and 
three different values of a, 0.9, 1.4 and 1.9 (slightly below, at, and slightly above the 
critical point, respectively). Right: Rescaled avalanche-size distribution showing good 
finite size scaling. This implies that the cut-off for the critical value (see Left) shifts 
progressively to the right, in a scale invariant way, upon enlarging the system size. 



to reproduce the value J^^'* = 0.025 for = 300 used in simulations in [29]. Results are 
mostly insensitive to the choice of P^^. 

Running computer simulations of Eq. (|3]) with this set of parameters, a steady state 
for both Vi and Jij is eventually reached, after an initial transient. In such a regime, 
driving events generate avalanches of activity. Fig. [1] shows time series in the steady 
state for the network-averaged value of uJij, uJ, with 

Results correspond to A^ = 1000 and three different values of a, 0.9, 1.4 and 1.9. 
Large avalanches (which are much more frequent in the supercritical phases) correspond 
to abrupt falls in uJ, while in between avalanches uJ grows linearly in time owing to 
the external driving. 

Observe the intermittent response of the system in all cases: peaks of activity of 
various sizes appear in all cases; note also the "quasi-periodic" behavior in all the three 
cases (similar quasi-periodic behavior had already been described for the Markram- 
Tsodyks model [321 [36] ). 

In order to determine the critical point, in Fig. |2]^left) we show the associated 
avalanche-size distributions for the same three values of a. All of them show, for small 
values of s a power-law decay, with exponent close to 1.5; for a = 0.9 (subcritical) there 
is an exponential cut-off while for a = 1.9 (supercritical) there is a "bump" for large 
size values, which defines a characteristic scale. In the intermediate case, a = 1.4 there 
is also an exponential cut-off but, upon increasing system size, it shifts progressively to 
the right in a scale invariant way, as corresponds to a critical point. This is illustrated 
in Fig. [2]^right) where critical distributions (i.e. for a = 1.4) for various system-sizes 
have been collapsed into a unique scale-invariant curve. 
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In Fig. m^top) we plot the distributions of uJ for different values of a, obtained by 
sampling values of J all along the dynamics. Observe the progressive broadening and 
displacement to the right upon increasing a. 
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Figure 3. Top: Probability distribution of uJ for a system size N — 1000 and different 
values of a. Only for a > ac — 1-4(1), the right tail of the distribution extends beyond 
the critical value of the static model af (A^ = 1000) uJc = 0.95. Bottom: P{uJ) 
at the critical point, ac — 1.4, for different system sizes; the width of the distribution 
does not decay with increasing system size and, therefore, this distribution is not delta- 
peaked in the thermodynamic limit. This reflects the fact that, for sufficiently large 
values of a the system hovers around the critical point alternating subcritical and 
supercritical regimes. For smaller values of a the system is always subcritical. 

Fig. [3]^bottom) illustrates the presence of strong finite size effects; in particular, for 
the critical point a = 1.4, we see that the distribution of J moves progressively to the 
right. The main observation to be made is that these distributions do not converge to 
narrower ones upon enlarging system size. Similar broad distributions are typical of non- 
conserving self-organized models, for which delta-peaked distributions are not obtained 
even if the infinite-size limit is taken This means that the dynamical model does 
not correspond to the static one with some fixed "effective" or averaged value of J, 
but to a dynamical convolution of different values of J, distributed in some interval 
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Figure 4. Phase diagram for the LHG model for different system sizes. Observe 
the presence of a critical line separating an active (supercritical) from an absorbing 
(subcritical) phase. Also, for large values of a a non-stationary or "explosive" phase 
(in which potentials grow unboundedly) exists. 



[Jmin{ci), Jmax{ci)], with Weights giveii by the distributions above. The probabihty of 
finding the system at any point out of such an interval [Jmm(tt), Jmax{c()] is zero (within 
numerical precision). 

As illustrated in Fig. H] we have verified that for u{J) > 1 (which occurs for 
large values of a; in particular, for a — > oo when N ^ oo limit), the loading 
mechanism dominates over the discharging one (synaptic depression), and the potential 
grows unboundedly with never ceasing activity; this is a non-stationary supercritical or 
explosive phase, analogous to the one reported for the static version of the model. 

Finally, we have also computed the average value of J at spike, i.e. right before 
the corresponding pre-synaptic neuron fires and before the value of J is diminished 
(see Fig. [5]). This quantity, that we call Jgp, appears in the analytical approach 
to be discussed below. Observe in Fig. [5], in analogy with the histograms above, 
the existence of broad distributions whose width does not decrease significantly upon 
enlarging system-size. Analyzing the highly non-trivial structure of these (multi-peaked) 
histograms is beyond the scope of this paper, but let us just mention that similar 
histograms with various peaks appear in related non-conserving model of SOqC |35j . 
Note also that they extend beyond uJsp = 1, even if their average is close to unity. 

3.2.2. Characterization of criticality Perusal of either Fig. [1] or Fig. [3t^top) leads to the 
important observation that it is only for values of a above the critical point (oc ~ 1.4) 
that the support of the distribution of uJ extends beyond the (A^-dependent) critical 
point of the static limit, i.e. that uJ^ax > af"*''^(A^) (see Fig. [3]^up)). For a < ac 
the dynamics is subcritical at every time (i.e. uJmax{N) is always below the threshold 
of the static model, q;^*"**'^(A^)), and hence avalanche distributions, being a dynamical 
convolution of avalanches with instantaneous subcritical parameters, are subcritical. 
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Figure 5. Probability distribution of values of uJ computed at spike, i.e. at their 
local maxima, just before being depressed. Curves correspond to different system sizes 
(from 300 to 1000) and fixed a, a — -i > Oc, i.e. in the supercritical phase. Observe the 
broad distribution, whose width does not decrease significantly upon enlarging system- 
size. Similar broad histograms, typical of SOqC systems, are obtained for other values 
of a. 



Instead, for a > uJ^ax > a^*"**'^(A^), and one can observe instantaneous values of 
the average synaptic strength, J, above the static model critical point, giving raise to 
instantaneous super-critical dynamics and system- wide propagation (observe that, in a 
fully connected topology, any site/neuron can be reached within one time-step). Then, 
during the avalanche, owing to the term —uJij5{t — t{^) in the second equation of Eq. ([3]), 
uJ decreases, and the system moves progressively from the supercritical regime to the 
subcritical one. This, in turn, becomes supercritical again upon recovering/loading. 
This cyclical shifting (analogous to the one for SOqC as described in [1]) provides a 
dynamical mechanism for the generation of a broad distribution of avalanche sizes in 
the steady state. 

Thus, it is only for a > Uc that arbitrarily large avalanches appear, and the critical 
point of the dynamical model corresponds, for any system size, to the value of a for 
which the maximum of the support of the distribution of values of J, Jmax, coincides 
with the critical point of the static model: 

Fig. [6](left) illustrates the coincidence (within numerical resolution) of the critical line 
for the static model, uJc{N), and the maximum of the support of the distribution of J at 
the critical point of the dynamical model for various system sizes. The small deviation 
between the two curves stems from the binning procedure employed to determine Jmax- 
The average value of J at the critical point is also plot in Fig.[6]^left) for illustration: 
at criticality, the average value is always far below unity, i.e. far below the conserving 
limit. Even in the infinite size limit, this curve remains below 1 (as a consequence of 
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Figure 6. Left: Critical value of J, Jc, in the static model (upper curve), maximum 
of the support of the distribution of J, Jmax, in the dynamic model (central curve), 
and average value of J, i.e. (J) at the critical point (lower curve). Note that this last 
curve lies in the subcritical region: (J) is not equal to 1 at the critical point. The 
dashed bell-shaped curve represents in a sketchy way the J-probability distribution for 
N — 2000; its height is unrelated to the a;-coordinate in the main graph; the peak is 
located around 0.88 (in coincidence with the (J) curve), while the upper tail of the 
distribution "touches" the vertical line, around 0.95 (i.e. at the corresponding point in 
the Jmax curve). Right: Scaling of the distance to the infinite-size critical point (i.e. 
1) in both the static and the dynamic LHG model as a function of the system size. As 
predicted by the general theory for non-conserving self-organized models, they both 
are power-laws with an exponent close to 1/3 (dashed line). 



the fact that the maximum of the distribution converges to 1 and the distribution is not 
a delta-function). 

Using our numerical estimates of the critical point as a function of (taken from 
Fig. m^left)) we have shown (see Fig. [6]^ right)) that the critical point converges to unity 
as 1 - uJ^axiN) ~ Ar-o-36(6). The same property holds for the static model, for which 
we obtain 1 — uJc{N)) ~ , This illustrates that the progressive shifting of 

the distributions in Fig. [3]^bottom) to the right occurs at the same pace as that of the 
critical point of the static model, in such a way that our estimate of the critical point, 
ac, is hardly sensitive to finite-size effects: for every studied system size, we obtain 
ac{N) 1.4(1) as illustrated in Fig. H 

Using the absorbing state picture of non-conserving self-organized systems, which 
predicts scaling to be controlled by a dynamical percolation critical point, we made the 
quantitative prediction that, for generic SOqC systems, the finite-size correction to the 
critical point should scale with system size as A^~^/^ (see and Section 5 below). In 
our case, 

uJrr.axiN ^ Oo) - uJ,nax{N) ~ N-'^^ (8) 

in agreement with the numerical estimates (see the dashed line in Fig. [6|^right)). This 
supports the validity of the theoretical framework presented in |1] to account for the 
present model: the critical behavior of neural avalanches is controlled by a dynamical 
percolation critical point. 
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Figure 7. Propagation of activity as a function of time in a two-dimensional (100* 100) 
implementation of the LHG model, for a = 2, in the supercritical phase. 



Before finishing this section, let us briefly present some results for the LHG model 
implemented on different type of topologies. In particular, we have numerically studied 
a version with a finite connectivity (random neighbors) as well as a two-dimensional 
lattice. In both cases, we find sub-critical and super-critical phases separated by a 
critical point, as in the fully connected lattice. 

For the random neighbor case, some details, as the way the cut-offs scale with 
system size, are different, but the main results are as in the fully connected case. 

For the two-dimensional lattice. Fig. [7] illustrates the evolution of an avalanche of 
activity for a particular value of a (in the supercritical regime). Observe the presence of 
a noisy wave of activity propagating outward from the seed; similar avalanches cannot be 
visualized in the fully connected case where activity reaches all sites in a single time-step. 
The waves shown in Fig. [7] resemble very much the ones observed in the retina (which is 
an almost two-dimensional network) before maturation [37] and, more importantly for 
the discussion here: they are fully analogous to supercritical waves appearing in 

• other non-conserving self-organized systems as forest-fires and 

• the dynamical percolation theory in the supercritical regime. 

This observation confirms, once again, the very close relationship between the LHG 
model and the theory of SOqC [1]. 

In summary, the LHG model is a representative of the class of non-conserving 
self-organized systems or SOqC, in which, as shown in a previous paper exhibits a 
conventional critical point separating a subcritical from a supercritical phase. Criticality 
is controlled by the maximum of the support distribution of J, Jmax, and not by 
its average value. This is in accordance with the general criterion for criticality in 
SOqC systems put forward in [1]: criticality emerges when the temporarily changing 
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background field, J, overlaps with the active phase of the underlying (static) absorbing 
state phase transition |3]. This result is of relevance for the analytical approach in the 
next section 



4. Analytical results 

The main conclusion of the previous section, i.e. the need to fine tune a to observe 
true criticality, seems to be in disagreement with the one presented in [29] for the LHG 
model. There it was claimed, relying on a mean-field calculation, that all values of a 
in the interval [l,oo[ are strictly critical. Using the hindsight gained from the results 
above, it is not difficult to find where the problem lies, as we show in what follows. 

Let us first construct (following LHG) a balance equation for the static limit of 
the model. Calling A*^* the inter-spike interval (time between two consecutive spikes 
of a given neuron) and A**^* the inter-avalanche interval (time between two consecutive 
avalanches, started at any neuron), the average number of avalanches between two spikes 
of the same neuron, (M), is 

(M) = (9) 

Obviously, A**^* has to be inversely proportional to P^^: if the external driving is reduced 
by a factor r the average time to generate an avalanche grows by a factor r. Levina et 
al. [29] actually showed that 

Jext ' ^ ' 

where, as above, e(A^) vanishes for — )■ oo. Focusing on a single neuron, in the steady 
state, it must obey the following balance equation 

Text J 

V^^^ = —A-^ + j^{s){M), (11) 

which equates the potential decrease for each spike (l.h.s. term) to the total potential 
increase between two consecutive spikes; this comes from two possible sources: (1) the 
average loading owing to external driving between two consecutive spikes (first term 
in the r.h.s.) and (2) the average charging from avalanches (second term). Note that 
— 1 is the number of neighbors of a given neuron and (s) is the averaged avalanche 
size. Fixing Vmax = 1, e(^) = 0, and plugging Eq. ([9]) and Eq. (ITOl) into Eq. (fTTI) . one 
readily obtains 

A^ 

ocmJ(s) (12) 



/\isi Jext 

for large values of A^. 

In the static case, {s) is given by Eq. ([5]), so A*^* can be expressed as a function of 
J, A^ and /^^'*. We have numerically verified that the resulting balance equation holds. 

On the other hand, in the dynamical case, J is not a constant and we do not have 
a simple expression for (s). The authors of [29] assume that the average avalanche 
size can still be written using Eq. but replacing uJ by u{Jsp). In particular, it is 
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hypothesized that avalanches can be effectively described as static avalanches with an 
effective branching rate given by the average branching ratio at spike (i.e. the synapses 
which are about to spike are the ones controlling the branching process of activity); this 



This equality is expected to hold in the infinite system-size limit and for infinitely 
large avalanches (where the law of large numbers applies) in which case, the average of 
sampled values of Jsp along sufficiently large avalanches can be safely replaced by (Jsp). 
In any case, it can be valid for only for branching ratios up to 1 (for which the geometric 
series converges). Substituting Eq. ( !T3l) into Eq. (fT2l) LHG readily obtain 



which, trivially, is smaller or, at most, equal to 1. From this, one concludes that the 
effective branching process is either subcritical or critical, but cannot be super-critical. 
Two comments are in order: 

The first one is that Eq. ( !T2|) is valid if and only if u{Jsp) is not larger than 1, hence, 
the calculation above does not exclude the existence of other (super-critical) solutions, 
for which Eq. f|T2l) would not hold. Actually, as illustrated in the numerics, for any 
finite system, an exploding phase, with branching ratio larger than unity, does exist (as 
a mater of fact, given a fixed value of a, depending on how the "loading" constant tj is 
scaled with system size, i.e. depending on how fast is the recovery of synapses, one can 
shift the location of the critical point and enlarge or reduce the size of the supercritical 
region) . 

The second one is as follows: the main approximation of the calculation above 
is the replacement of the average of sampled values of Jsp along any sufficiently large 
avalanche by {Jsp)- If, during avalanche propagation, the uncovering of values of Jsp 
from P{Jsp) (which is depicted in Fig.[5]for a particular value of a) occurred in a random, 
uncorrelated, way then the process would be what is called in the literature a "branching 
process in a random environment" [38]. Such a process turns out to be controlled by 
the average value of the random branching ratio ^8]. In such a case, the calculation 
would be exact and, for any value of a for which the average branching ratio is unity, 
the process would be critical. 

However, the uncovering of values of Jsp in the LHG model exhibits strong 
correlations. Jsp fluctuates around the central value uJsp = 1 in a rather correlated way. 
This is illustrated in Fig. [HI where we plot a return map for of uJsp averaged along each 
single avalanche. Notice that the return map is not structureless as would correspond 
to a random process. Instead, the system is progressively charged towards large values 
of the synaptic intensity and, afterwards, it gets suddenly discharged, starting a new 
cycle. In this way, the true dynamics of the system consists of a continuous alternation 
of supercritical (where most of the Jsp take values above 1), and subcritical dynamics: 
individual avalanches are either subcritical (average branching ratio smaller than 1) or 



is: 




(13) 




]Y2 _ ^isijext 



(14) 
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Figure 8. Return map for uJgp averaged along two consecutive avalanclies An and 
An+i in the supercritical regime. The broken line joints (clockwise) 20 consecutive 
points of the map to illustrate the temporal structure of the charging-discharging cycle. 
The non-trivial structure of the map reflects the presence of strong correlations: the 
system typically moves up in a few steps along the main diagonal (see the broken line) 
then, after reaching the supercritical regime uJsp > 1, a large avalanche is produced, 
and the system returns back to the subcritical regime uJgp < 1, to start a new charging- 
discharging cycle. The diagonal dashed-line, uJsp{An+i) — uJsp{An), is plotted as a 
guide to the eye. 

supercritical (average branching ratio above 1), and hence the resulting avalanche-size 
distribution is a complex one (not a simple power- law). This is illustrated in Fig. ([9]) 
which shows the avalanche size distribution for different system sizes and a = 4 which lies 
in the supercritical phase (the rest of parameters are as in Fig. |5]). Even if the averages of 
uJsp (as calculated from Fig. [5]) are very close to 1 for all sizes, the curves in Fig. |5]show a 
bump at large avalanche-sizes, reflecting the presence of many supercritical avalanches. 
This effect does not decrease upon increasing system size, even if the bump moves 
progressively to larger values as the system size is increased. Similarly, correlations are 
also responsible for the shift from the predicted mean-field critical point a = 1 to the 
actual one ac ~ 1.4. 

In conclusion: even if the branching ratio turns out to be always very close to unity 
for any value of a > ac, the avalanche-size distributions are not generically pure power- 
laws. In the supercritical phase there are bumps revealing its non-truly scale-invariant 
nature. This conclusion is in agreement with the general scenario for non- conserving 
self-organized systems introduced in [1]. 

As a final remark, we want to emphasize that the results by Levina et al. are mostly 
correct: the branching ratio is actually equal to unity in a broad region of parameter 
space (in the infinitely large system size limit). However, as explained above, this 
"marginality" of the averaged branching ratio does not exactly corresponds generically 
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Figure 9. Avalanche size distribution for a = 4 (rest of parameters, as in plots above) 
and different system sizes (as in Fig. [S|) Observe the presence of bumps, which do not 
disappear by increasing system size. This illustrates the existence of a supercritical 
phase in the LHG model. 

to true scale invariance. 

The main virtue of the LHG model is that, even if not generically critical, it 
generates a rather broad "pseudo-critical" region, exhibiting partial power-laws. The 
ultimate reason for this is rooted in the extremely slow loading (recovering) process 
of the background field (synaptic strength), which occurs during avalanches. This is 
to be compared with the more abrupt loading in forest-fire and earthquake models, 
which occurs between avalanches. This more abrupt loading induces excursions around 
the critical point to be broader than their counterparts in the LHG model. This is 
particularly true in the (probably unrealistic) case in which tj diverges with system 
size. 

5. A simple absorbing-state Langevin equation approach 

In order to have a more explicit connection between the LHG model and the family 
of SOqC models and theory discussed in [4j, in this Section we construct a Langevin 
equation for the LHG model, which includes absorbing states (and is therefore a natural 
extension of the Langevin theory for SOC, as introduced in [21 [39l HQ]) and turns out 
to be almost identical to the general Langevin theory for SOqC systems (as introduced 
ini). 

For the sake of simplicity, and without loss of generality, let us consider 
homogeneous initial conditions for all Vi and Jij, i.e. Vi = V Wi and Jjj = J Wi,j, 
and also, /^^* = 0. Under these conditions, and given the deterministic character of the 
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dynamics, all neurons evolve synchronously and Eqs. (E]) can be simply rewritten as: 




Tj \U 

where tsp are the firing times. Let us remark that in order to treat the more general 
heterogenous case it suffices to keep sub-indexes in the different variables. 

The spike terms, proportional to 6{t — tgp), can be alternatively written as 

6it-tsp)^ P = QiV-V,na.), (16) 

where 6(x) is the Heaviside step function (we take the convention 6(0) = 0); i.e. spike 
terms operate only whenever the potential is above threshold, implying that the activity 
variable, p, is non-zero only in such a case. Thus: 



, (17) 
- uJp. 

Further analytical progress can be achieved by regularizing the step-function in Eq. flTBl) 
as a hyperbolic-tangent: 

p^^{l + tanh[/3{V -V^a.)]), (18) 




which is a good approximation provided that /3 ^ 1. Inverting Eq. f[T8|) : 

_ arctanh (2p - 1) + Vmax ,^ 
V - ^ (19) 

and, taking derivatives on both sides, 

which is well-defined provided p g]0, 1[. Let us underline that the forthcoming equations 
are also valid at p = 0, where activity ceases. Using this, Eq. f|T7|) can be rewritten as: 

dtp = 2/3 {uJ - Vmax) P^ (1 - p) 

dtJ= ^(^-j)-uJp ^''^ 



Tj \ U 

which, omitting higher order terms, reduces to: 



dtp= 2l3uJp'-2l3VmaxP' 

Tj\U 
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Renaming variables as: 2/3Vmax b, 2(3u — )■ tu, J — )■ 0, — ^7, Cij/u — )■ 0c, and 
u W2, one obtains: 

dtp = w(t)p^ - hp^ 
dt(f) = {(pc - 4>) - W2(f)p. 

The equation for p is a typical mean-field equation for a system with absorbing states 
(i.e. all dynamics ceases when p = 0). It includes a coupling term with the background 
field 0: the larger the background the more activity is created. In the simplest possible 
theory of SOqC (see [1]), such a coupling is linear in p, but the effect of both types 
of coupling can be argued to be qualitatively identical. On the other hand, the second 
equation is identical to the mean-field background equation for SOqC systems: the 
presence of activity reduces the background field while the loading mechanism, acting 
independently of activity, increases it. 

Except for the coupling term which is quadratic in p, these mean-field equations are 
identical to the ones proposed in |1] to describe non- conserving self-organized models at 
a mean-field level. Also, in analogy with SOqC systems, when slow driving is switched 
on, i.e. 7^ 0, activity can be spontaneously created, even if p = 0, generating 

avalanches of activity. Moreover, if some sort of stochasticity (and hence heterogeneity) 
is introduced into the dynamics, then it can be easily seen that: 

• a noise term, proportional to ^/p, needs to be added to the first equation, 

• a diffusion term accounting for the coupling with nearest neighbors, and 

• a linear-coupling term is perturbatively generated in the activity equation (and, 
thus, the quadratic-coupling becomes a higher order term). 

Therefore, after including fluctuations and omitting higher order terms, the final 
set of stochastic equations that we have derived is identical to the one of dynamical 
percolation [H] in the presence of a "loading mechanism" , i.e. to that of SOqC systems 
as described in 

This heuristic mapping between the LHG model and the general theory of SOqC, 
justifies from an analytical viewpoint all the findings in previous sections (including the 
quantitative prediction for the finite-size scaling of (1 — uJmax{N))) and firmly places 
the LHG model in the class of self-organized quasi-critical models, lacking true generic 
scale-invariance . 

6. Conclusions 

Cortical avalanches, first observed by Beggs and Plenz [5] , were claimed to be generically 
power-law distributed and, thus, critical. Such a claim led to an outburst of activity in 
Neuroscience trying to understand the origin and consequences of such a generic scale- 
invariance. At a theoretical level, Levina, Herrmann, and Geisel [29] proposed a simple 
model (a variation of the Markram-Tsodyks model for chemical synapses), claimed to 
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reproduce generically scale-invariance. In particular, these authors performed a mean- 
field calculation leading to the conclusion that, for any value of the control parameter, 
a, larger than unity, generic critical behavior is observed. They also conducted some 
computational studies to support their findings. 

The LHG model turns out to be very similar to slowly driven models of self- 
organized criticality such as earthquake and forest-fire models. As in these other models, 
and in contrast to sandpiles, in the LHG one the dynamics is non-conserving (refiecting 
the leaking/dissipative dynamics of actual synaptic signal transmission). 

It is by now a well-established fact that non- conserving self-organized models are 
not generically critical but just "hover around" the critical point of an underlying 
absorbing phase transition, with finite excursions (of tunable amplitude) into the active 
and the absorbing phases. As they do not converge to the critical point itself, generic 
scale-invariance cannot be invoked (see ^ and references therein). The term self- 
organized quasi criticality (SOqC) has been proposed to refer to such a class of systems, 
emphasizing the differences with conserving SOC models. 

Given the contradiction between this general result and the claim in [29], in this 
paper we have scrutinized the LHG model, both numerically and analytically, and 
reached the following conclusions: 

• Both in its static and its dynamical form, the model exhibits absorbing and active 
phases and a non-trivial critical point separating both of them. 

• It is only if parameters are fine tuned to such a critical point that true scale- 
invariance emerges and the distribution of avalanche-sizes is power-law distributed. 

• The mean-field calculation in p9], supporting generic criticality, lead indeed to a 
branching ratio equal to unity in a broad interval of phase space, but this does not 
imply generic scale-invariance. 

• A Langevin equation, including absorbing states, has been derived for the LHG 
model. Such an equation reduces to the analogous one proposed to describe 
generically non- conserving self-organized (SOqC) models. Thus, all the general 
conclusions obtained from such a theory in [1] apply to the LHG model, providing 
analytical support to the numerical findings above. 

It is worth stressing that our results do not subtract merit from the LHG model. 
Actually, strict criticality might not be required to explain the truncated power-laws 
reported by Beggs and Plenz; the dynamical LHG model generates partial power- 
laws compatible with the empirical findings by Beggs and Plenz for a relatively broad 
parameter (a) interval, as shown in Fig. [TOl Moreover, the fact that the model can 
generate critical, subcritical, and supercritical regimes, depending on parameter values, 
converts the LHG model into an adequate one to describe the state-of-the-art in neuronal 
avalanches. As mentioned in the Introduction, Pasquale et al. have shown in a recent 
paper [28] that, depending on several experimental features, cortical avalanches can 
indeed be either critical, subcritical, or supercritical. 
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Figure 10. Range of compatibility between the results of the LHG model, for different 
values of N, and the empirical results by Beggs and Plenz; for large system sizes 
{N > 700) values of a between 1 and 1.4 give avalanche-size distributions compatible 
with those observed by Beggs and Plenz [5], even if they are subcritical. 



The main implication of our work can be summarized as follows: if future 
experimental research conducted on cortical networks were to support that critical 
avalanches are the norm and not the exception, then, one should look for more elaborate 
theories, beyond simple self-organization, to explain this. Standard self-organization 
does not suffice to explain criticality in non-conserving systems. Parameters have to 
be tuned or "selected" to achieve a close-to-criticality regime. For instance, the claim 
by Royer and Pare [31] that homeostatic regulation mechanisms keep cortical neural 
networks with an approximately constant (i.e. conserved) global synaptic strength could 
be at the basis of such a less generic theory beyond simple self-organization. Another 
inspiring possibility is that natural selection by means of evolutionary and adaptive 
processes leads to parameter selection, favoring critical or close-to-critical propagation 
of information in the cortex [30] . A more realistic approach should also include long-term 
plasticity [12] , as well as co-evolutionary mechanisms, shaping the network topology. We 
shall explore these possibilities in a future work. 

Appendix : Synchronization and oscillatory properties 

Synchronization was studied in the self-organized criticality literature as a possible 
mechanism, alternative to conserving dynamics, leading to generic scale invariance |13] . 
Even though such a suggestion turned out not to be true [1], let us explore here the 
oscillatory and synchronization properties observed in numerical simulations of the LHG 
model. With this aim, we compute the power-spectra, S{f), for the time series of J 
shown in Fig. [1] (as well as for other values of a). In all illustrated in Fig. [TTl 

the spectra exhibit peaks at some characteristic frequencies, /. Closer inspection reveals 
that the maximum peak appears at a characteristic frequency, /c, which we have verified 
to be inversely proportional to (A*'^*). This indicates that the typical time needed for a 
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Figure 11. Power spectrum of the LHG model for a = 2 (i.e. supercritical). 
Frequencies / are plotted rescaled by a factor (A***). Note the presence of peaks, at 
/ cx (A***)~^, coexisting with fat tails. The tail decay fc~^ (dashed line) is characteristic 
of sawtooth profiles (i.e. with linear increases). 



neuron to overcome threshold and spike again introduces a characteristic scale into the 
system, entailing periodicity. Observe also that the power-spectra exhibit fat tails, with 
exponent k~^, characteristic of sawtooth profiles with linear increases. 

The previous numerical analysis can also be done for the static model, with almost 
identical results: the origin of the periodic behavior lies in the charging/discharging 
cycle of potentials, V, and is not crucially affected by the synaptic strengths being fixed 
or not. 

Given that individual neurons oscillate with a certain periodicity, let us 
study (in analogy with other analyses of non- conserving self-organized systems) the 
synchronization (or absence of it) between different units (either neurons or synapses). 

In order to quantify synchronization, we use bins of size 2 ■ 10"'' (for V) and 10~^ 
(for J), and consider as an order parameter the fraction of neurons/synapses which are 
synchronized, i.e. which lie in the same discrete bin, divided by the number of occupied 
bins. Such a parameter becomes arbitrarily small for a large enough random system and 
is 1 in the case of perfect synchronization. If the total number of elements into a multiply 
occupied bin is Ns and the number of bins is Nb, the value of the synchronization order 
parameter, 0, is 

<Pv-^ 0, (24) 

for neurons {V) and synapses (J), respectively. By monitoring (pv, we observe that 
the potentials in the system converge to a totally un- synchronized state (0y ~ 10~^). 
This is in agreement with the uniform distribution of values of V employed in analytical 
arguments above. 

On the other hand, by measuring 0j we observe that it rapidly converges to a 
stationary state value, Ni, = N, reflecting a perfect synchronization of the different 
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synapses of any given neuron, j, i.e. Jij = Jkj for any values of i and k: all the synapses 
Jij emerging from of a given (presynaptic) neuron, j, converge to a common state. This 
can be easily understood using the following argument. The dynamics of Jij and Jkj 
are controlled by the same equation 



dJ, 



1,3 



-(--J,,)-M,,r/(t), (25) 



dt Tj 

where / is either i or k and ri{t) is a (positive) noise, accounting for the spikes of (pre- 
synaptic) neuron j, which is obviously common to all synapses of j. Subtracting Eq. (!25|l 
for k from Eq. ( !25ll for i we obtain that the difference, A = Jij — Jkj evolves as 

^ = -A[- + nr/(t)], (26) 

which, given the positivity of tj, u and r], entails a negative Lyapunov exponent and, 
hence, convergence to the synchronous state, Jjj = J^j - i-e. all synapses emerging 
from a given pre-synaptic neuron synchronize. Observe that this type of synchronization 
is similar (but not identical) to that observed in, for instance, earthquake models 
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